#include "image.h"

image::image() {
 
  ifstream myimage;
  vector<double> temp;
  string input;
  int t;
  double s;
  epsilon = 0.0001; //arbitrary

  cout << "Please enter input file: " << endl;
  cin >> input;
  myimage.open(input.c_str());
  
  if (myimage.is_open()) {
    while(myimage.good()) {
      for (int i=0;i<2;i++) {
	myimage >> t;
	size.push_back(t);	
      }
      temp.clear();
      for (int i=2;i<size[0]+2;i++) {
	temp.clear();
	for (int j=0;j<size[1];j++) {
	  myimage >> s;
	  temp.push_back(s);
	}
	initial.push_back(temp);
      }
    }
    myimage.close();
  }
  else {
    cout << "Unable to open image." << endl;
  }
  
}

image::~image() {

}

void image::print(vector < vector <double> > matrix, int m, int n) {
  if (!matrix.empty()) {
    for(int i=0;i<m;i++) {
      for(int j=0;j<n;j++) {
	cout << matrix[i][j] << "\t";
      }
      cout << endl;
    }
  }
  else {
    cout << "No file has been read." << endl;
  }
}

void image::print3D(int size) {

  if(!ThreeDmatrix.empty()) {
    for(int i=0;i<size;i++) {
      for(int j=0;j<size;j++) {
	cout << ThreeDmatrix[i][0][j] << "\t";
      }
      cout << endl;
    }
  }
}

vector <vector <double> > image::computetranspose(vector < vector <double> > matrix, int m, int n) {

  vector <double> temp;
  vector < vector <double> > tempmatrix;
  double t;
  temp.clear();
  tempmatrix.clear();
  
  for(int i=0;i<n;i++) {
    for(int j=0;j<m;j++) {
      t = matrix[j][i];
      temp.push_back(t);
    }
    tempmatrix.push_back(temp);
    temp.clear();
  }

  return tempmatrix;

}

vector < vector <double> > image::multiply(vector < vector <double> > matrix1, vector < vector <double> > matrix2, int m, int n, int p, int q) { //multiplies matrix 1 (mxn) by matrix 2 (pxq)

  vector <double> tempvector;
  vector < vector <double> > tempmatrix;
  double t;
  tempvector.clear();
  tempmatrix.clear();

  if(n==p) {
    for(int i=0;i<m;i++) {
      tempvector.clear();
      for(int j=0;j<q;j++) {
	t = 0;
	for(int k=0;k<n;k++) {
	  t = t + matrix1[i][k]*matrix2[k][j];
	}
	tempvector.push_back(t);
      }     
      tempmatrix.push_back(tempvector);
    }
  }
  else {
    cout << "Cannot multiply matrices, incompatible size." << endl;
  }
  return tempmatrix;
}

vector < vector <double> > image::subtract(vector < vector <double> > matrix1, vector < vector <double> > matrix2, int m, int n, int p, int q) { //matrix 1 minus matrix 2

  vector <double> tempvector;
  vector < vector <double> > tempmatrix;
  double t;
  tempvector.clear();
  tempmatrix.clear();

  if(m == p && n ==q) {
    for(int i=0;i<m;i++) {
      tempvector.clear();
      for(int j=0;j<n;j++) {
	t = 0;
	t = matrix1[i][j]-matrix2[i][j];
	tempvector.push_back(t);
      }
      tempmatrix.push_back(tempvector);
    }
  }
  else {
    cout << "Cannot subtract matrices, incompatible size." << endl;
  }
  return tempmatrix;
}

bool image::check_tolerance(vector < vector <double> > matrix, int m, int n) { //checks to see if each entry is within epsilon of 0

  for(int i=0;i<m;i++) {
    for(int j=0;j<n;j++) {
      if (matrix[i][j] > epsilon || matrix[i][j] < epsilon*(-1)) {
	return false;
      }
    }
  }
  return true;
}

void image::computemodulus() { //computes modulus THIS DOESNT WORK YET

  vector <double> temp;
  double t;
  temp.clear();
  modulus.clear();

  for(int i=0;i<size[1];i++) {
    temp.clear();
    for(int j=0;j<size[1];j++) {
      t = sqrt(hermetian[i][j]);
      temp.push_back(t);
    }
    modulus.push_back(temp);
  }
  return;
}

double image::dotProduct(vector < vector <double> > v1, vector < vector <double> > v2, int size) { //only for use in eigenvalue computation function

  double temp = 0;

  for(int i=0;i<size;i++) {
    temp = temp + v1[i][0]*v2[i][0];
  }

  return temp;

}

vector < vector <double> > image::normalize(vector < vector <double> > v, int n) {

  double size = 0;
  vector < vector <double> > tempmatrix;
  vector <double> tempvector;
  double t = 0;
  tempmatrix.clear();
  tempvector.clear();
  
  for(int i=0;i<n;i++) {
    size = size + v[i][0]*v[i][0];
  }
  size = sqrt(size);
  
  for(int i=0;i<n;i++) {
    t = 0;
    t = v[i][0] / size;
    tempvector.push_back(t);
  }

  tempmatrix.push_back(tempvector);  
  tempmatrix = computetranspose(tempmatrix,1,n);
  return tempmatrix;

}

void image::create3DArray(int size) {

  double temp = 0;
  vector <double> tempvector;
  vector < vector <double> > tempmatrix;

  srand(unsigned(time(NULL)));

  for(int i=0;i<size;i++) {
    for(int j=0;j<size;j++) {
      temp = ((double)rand()/(double)RAND_MAX);
      tempvector.push_back(temp);
    }
    tempmatrix.push_back(tempvector);
    tempvector.clear();
    ThreeDmatrix.push_back(tempmatrix);
    tempmatrix.clear();
  }

}

vector < vector <double> >  image::gramSchmidt(int size) {

  vector < vector <double> > tempmatrix;
  vector <double> tempvector;
  double tempvalue;
  double tempvalue2;
  vector < vector < vector <double> > > temp3Darray;
  tempmatrix.clear();
  tempvector.clear();
  tempvalue = 0;
  tempvalue2 = 0;

  for(int i=0;i<size;i++) {   
    tempmatrix.clear();
    tempvector.clear();
    tempvalue = 0; 
    for(int j=0;j<size;j++) {
      tempvalue = ThreeDmatrix[i][0][j]; //check for permutation?
      tempvector.push_back(tempvalue);
    }
    tempmatrix.push_back(tempvector);   
    tempvalue = 0;
    for(int j=0;j<i;j++) {
      tempvector.clear();
      tempvalue = 0;
      tempvalue2 = 0;
      for(int k=0;k<size;k++) {
	tempvalue = tempvalue + tempmatrix[0][k]*ThreeDmatrix[j][0][k];
	tempvalue2 = tempvalue2 + ThreeDmatrix[j][0][k]*ThreeDmatrix[j][0][k];
      }
      tempvalue = tempvalue / tempvalue2;
      for(int m=0;m<size;m++) {
	ThreeDmatrix[i][0][m] = ThreeDmatrix[i][0][m] - tempvalue*ThreeDmatrix[j][0][m];
      }
    }
  }

  //normalize
  for(int i=0;i<size;i++) {
    tempvalue = 0;
    for(int j=0;j<size;j++) {
      tempvalue = tempvalue + ThreeDmatrix[i][0][j]*ThreeDmatrix[i][0][j];
    }
    tempvalue = sqrt(tempvalue);
    for(int j=0;j<size;j++) {
      ThreeDmatrix[i][0][j] = ThreeDmatrix[i][0][j] / tempvalue;
    }
  }

  //transpose
  tempvector.clear();
  tempmatrix.clear();
  tempvalue = 0; 
  for(int i=0;i<size;i++) {
    for(int j=0;j<size;j++) {
      tempvalue = ThreeDmatrix[j][0][i];
      tempvector.push_back(tempvalue);
    }
    tempmatrix.push_back(tempvector);
    temp3Darray.push_back(tempmatrix);
    tempmatrix.clear();
    tempvector.clear();
  }
  for(int i=0;i<size;i++) {
    for(int j=0;j<size;j++) {
      ThreeDmatrix[i][0][j] = temp3Darray[i][0][j];
    }
  }


  

}


vector < vector <double> > image::computeeigenvector(vector < vector <double> > matrix, int n) {

  eigenvector.clear();
  for (int i=0;i<n;i++) {
    eigenvector.push_back(1); //push in random value instead?
  }

  double t;
  vector <double> tempvector;
  vector < vector <double> > tempmatrix1;
  vector < vector <double> > tempmatrix2;
  tempmatrix1.clear();
  tempmatrix2.clear();
  tempvector.clear();
  t=0;

  tempmatrix1.push_back(eigenvector);
  tempmatrix1 = computetranspose(tempmatrix1,1,n);
  tempmatrix1 = normalize(tempmatrix1,n);
  tempmatrix2 = multiply(matrix,tempmatrix1,n,n,n,1);
  tempmatrix2 = normalize(tempmatrix2,n);
  
  while(!check_tolerance(subtract(tempmatrix2,tempmatrix1,n,1,n,1),n,1) ) {
    tempmatrix1 = multiply(matrix,tempmatrix1,n,n,n,1);
    tempmatrix1 = normalize(tempmatrix1,n);
    tempmatrix2 = multiply(matrix,tempmatrix2,n,n,n,1);
    tempmatrix2 = normalize(tempmatrix2,n);
  }
  
  print(tempmatrix2,n,1);
  return tempmatrix2;
}

double image::computeeigenvalue(vector < vector <double> > matrix, int n) {

  vector < vector <double> > tempmatrix;
  double eigenvalue;

  tempmatrix = multiply(matrix,testing,n,n,n,1);
  for(int i=0;i<n;i++) {
    if(testing[i][0]!=0) {
      eigenvalue = tempmatrix[i][0] / testing[i][0];
    }
  }

  return eigenvalue;
}
